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ABSTRACT 

This paper contains a local linear stability analysis for accretion disks un- 
der the influence of a global radial entropy gradient (3 = —dlogT/dlogr for 
constant surface density. Numerical simulations suggested the existence of an 
instability in two- and three-dimensional models of the solar nebula. The present 
paper tries to clarify, quantify, and explain such a global baroclinic instability 
for two-dimensional fiat accretion disk models. As a result linear theory predicts 
a transient linear instability that will amplify perturbations only for a limited 
time or up to a certain finite amplification. This can be understood as a result 
of the growth time of the instability being longer than the shear time which 
destroys the modes which are able to grow. So only non-linear effects can lead 
to a relevant amplification. Nevertheless, a lower limit on the entropy gradient 
oc /5 ~ 0.22 for the transient linear instability is derived, which can be tested in 
future non-linear simulations. This would help to explain the observed instability 
in numerical simulations as an ultimate result of the transient linear instability, 
i.e. the Global Baroclinic Instability. 

Subject headings: accretion, accretion disks — circumstellar matter — hydro- 
dynamics — instabilities — turbulence — methods: numerical — solar system: 
formation — planetary systems 
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1. 



Introduction 



In Klahr & Bodenheimer 2003 (also Klahr & Bodenheimer 2000) we presented two- 
dimensional (and three-dimensional) numerical hydrodynamical simulations of accretion 
disks, without magnetic fields and self-gravity. In the case with a radial entropy {S) gradient, 
the disk developed an instability, while in the case of no radial entropy gradient it stayed 
perfectly quiet. 



with c^:specific heat, P: vertically integrated pressure, S: vertically integrated density and 
7: the adiabatic index (P ~ E'''). It is necessary to define a 7 for the surface density, because 
the entire paper will use the two-dimensional fiat approximation of hydrodynamics, even it 
is just a crude approximation of the much more complicated relation between temperature 
and surface density. Goldreich, Goodman & Narayan (1986) derive 7 = (3r — l)/(r -|- 1) 
for this case, which relates to 7 = 1.354 for the three-dimensional adiabatic index F = 1.43 
which is realistic for the expected He, H2 mixture of the solar nebula. 

The non-axisymmetric instability in Klahr & Bodenheimer (2003) formed waves and 
grew exponentially at a rate of ~ (O.lJl) { Q — orbital frequency), eventually breaking up 
into stable vortices. The formation of vortices is of special interest, because vortices can be 
the trigger for a fast and very efficient planet formation process (Klahr 2003). 

Measurements of the Reynolds stresses indicated angular momentum transport radially 
outward at an efficiency comparable to a = 10~^ — 10^^ {a follows from the Shakura-Sunyaev 
definition a = with u efi^ective viscosity, Cg speed of sound and H the pressure scale 

Cs ^ 

height). In order to understand this instability a linear stability analysis is performed in this 
paper. The result is a critical value for the entropy gradient, and the growth rates as they 
depend on the wave number and time. 



Barotropic (P = P(S)) Keplerian disks are centrifugally stable as predicted by the 
Rayleigh criterion (see e.g Balbus, Hawley, & Stone 1996 for numerical and analytical in- 
vestigations). The situation changes if self gravity (Toomre 1964) or magnetic fields become 
important (Balbus & Hawley 1991). But both mechanisms are restricted to parameter ranges 
unreahstic for a large portion of protoplanetary disks (Gammie 1996). Therefore, one still 
looks for alternative effects that could drive turbulence and thus angular momentum trans- 
port in accretion disks. Riidiger, Arlt & Shalybkov (2002) showed that in vertically stratified 
disks there exists no exponentially growing linear mode without vertical shear, which was 
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suggested earlier. The baroclinic instability that is analyzed here in this paper does not 
consider the vertical structure explicitly, but only considers a radial entropy gradient. This 
might confuse readers who arc familiar with the Baroclinic Instability as dealt with in geo- 
physical contexts (e.g. Pcdlosky 1987). There, one considers the vertical structure to be 
important. But this is not entirely true, because one uses a variable vertical variation in 
density and velocity for the purpose of introducing a radial pressure gradient in a fluid that 
is conveniently assumed to be isothermal and incompressible otherwise. A linear analysis 
like the 2-layer model (see Pedlosky 1987) does not consider vertical modes. Thus, from now 
on we want to distinguish between the Classical Baroclinic Instability in the earth atmo- 
sphere, which formally arises if the potential vorticity gradient changes sign in the vertically 
stratified atmosphere, and the Global Baroclinic Instability considered here. What both 
instabilities have in common is the radial gradient of entropy, which makes the equilibrium 
state baroclinic (i.e. there is an inclination between pressure and density gradient). For both 
instabilities the radial entropy gradient is the sine qua non condition for instability. The 
main difference between the two instabilities is that Keplerian shear is not present in the 
Classical Barochnic Instabihty. Also, the 2-layer model for the Classical Baroclinic Insta- 
bility allows one to determine an upper cut-off for the azimuthal wavenumbers which are 
unstable; this is possible because the vertical shear kills the instability by eventually bringing 
the upper and the lower layers out of phase. This stabilizing effect of the vertical shear in 
accretion disks is not considered in this paper and will be investigated later. Tentatively, 
one can speculate that the vertical shear will damp high azimuthal wave numbers in disks. 
Those high wave number modes are also the strongest damped by the radial shear, as one 
will see in this paper. However, the paper focuses on the small to medium range of wave 
numbers anyway. 

The baroclinic instability in this paper can be regarded as an essentially convcctive 
instability in the radial direction but modified by rotation and shear. If there would be no 
shear but only rotation then axisymmetric perturbations obey the Solberg-H0iland stability 
criterion. But non-axisymmetrically, since the angular momenta of individual fluid elements 
is not conserved, the criterion reduces to the Schwarzschild-Ledoux criterion 4^4^ > as 

' ar clr 

shown by Cowling(1951). This criterion says that an instability occurs when radial pressure 
P and entropy S gradient point in the same direction. In some sense this paper is an 
extension of Cowlings work to apply it for accretion disks by incorporating the treatment of 
shear. 
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1.2. Outline of the paper 



Section 2 gives the basic hydrodynamic equations relevant for a two-dimensional flat 
accretion disk with a radial entropy gradient. In Section 3, a local linear analysis in La- 
grangian shear coordinates is performed. The result is a 4th order dispersion relation, which 
is a function of the time-dependent radial wave number kx{t) and the azimuthal wave number 
ky in local x — y coordinates. Section 4 expands the linear analysis to a global model, but 
the analysis is restricted to the radial modes kr{t = 0) = 0, showing that the existence of the 
non-axisymmetric instability depends on a minimum entropy gradient. A complete global 
analysis, as well as the determination of the eigenf unctions, is left to future work. Section 
5 contains a numerical solution of the linearized equations. Using the linearized equations 
from Section 3, it is possible to give an easy explanation of the special mechanism behind the 
Global Baroclinic Instability. This explanation is given in Section 6. Section 7 summarizes 
and discusses the results. 



Consider a two-dimensional flat disk in cylindrical coordinates. The vertically {z) inte- 
grated density is: 



The vertical velocity component (vz) vanishes at infinity, which means that there is no 
mass lost or gained from outside. The governing ideal hydrodynamic equations in polar 
coordinates are: 



2. Hydrodynamic Equations 




(2) 



(3) 



dtVr + VrdrVr H dfhVr 

r 



(4) 




(5) 




where P is the vertically integrated pressure. 




(7) 



and ^kep is the local Keplerian frequency. The energy equation. 
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is replaced by an equation for the pressure, Eth — P/{'y — 1), via the ideal gas equation with 
adiabatic index 7. The baroclinicity of a general flow is given by the baroclinic term: 

Vp(r, z, (p) X Vp(r, z, (p) ^ 0. (9) 

In a three-dimensional but axisymmetric equilibrium state of a disk, there is in fact a non- 
zero baroclinic term above and below the midplane, which vanishes at the midplane itself 
(z = 0): 

Vp(r, z) X Vp(r, z) ^ 0. (10) 

In the two-dimensional flat approximation used in this paper, this term vanishes for the 
axisymmetric mean state, 

VE(r) X VP(r) = 0, (11) 

but still the disk is baroclinic. This is shown by the existence of the two independent 
equations for density and internal energy, reflecting the inclination of the equipressure surface 
towards the cquidensity surface with radius. Again, most important for the instability is that 
the non-axisymmctric deviations from the mean state can lead to the rise of the baroclinic 
term even in two dimensions, 

VE(r,0)xVP(r,0)^O, (12) 

and vorticity can be generated. 



2.1. Defining the mean state 

The quantities Q = {T^.Vr.v^, P) can be split into a time independent axisymmetric 
mean part Q and a perturbation Q' as 

Q = Q{r) + Qir,<P,ty (13) 

The constant density background is assumed to follow a simple power law: 

SW = So(^) '\ (14) 
The pressure is distributed in the same fashion: 

^W = Po(^) \ (15) 
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thus for (5s the entropy gradient (Eq. 1) it follows: 

d^S = c,^^^^ and Ps = l3- (16) 
r 

There is no radial mean flow, which means: 

Vr = 0. (17) 

The mean rotation is axisymmetric and determined by the radial component of the Euler 
equation, thus the equilibrium between radial pressure gradient, gravity and centrifugal 
acceleration is: 

= -i9.P+^-QL.r, (18) 
which leads to the definition for the local orbital frequency: 



^-T-i^le^+^rP. (19) 



Using the definitions from Eqs.l4 and 15 one gets: 



r 

In the case of no radial pressure gradient 



/3E-/3-2 



^ — ^feep- (21) 

Unfortunately, Q does not in general follow a power law; but, as ^2^^^ ^ r"^, it is possible 
for /3s — /3 — 2 = —3 to define an Qq 

n^no(-) \ (22) 



To 



with 



^0- Xl^lep- -^kepxl"^- (23) 



This situation is for instance given if the surface density is constant over radius {jSj^ = 0) 
and the sound speed is a constant fraction of the Kcplcrian velocity; e.g., the relative 
pressure scale height H/r is constant with radius. This leads to the definition of two cases, 
in which A: there is no radial pressure gradient (/3 = 0) or B: a pressure gradient exists with 
P — 1. For all cases of our analytic treatment, we simphfy the density structure to he P-£ — 
(E = Eq) . In case A the rotational frequency is exactly Keplerian, but in case B the disk is 
slightly sub-Keplerian while still maintaining a Keplerian (~ r~^'^) profile. 
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3. Local Linecir Analysis 

This linear analysis will demonstrate the dependence of the instability on the exis- 
tence of non-axisymmetric modes and the radial entropy gradient. In contrast to the WKB- 
approximation for tightly wound spirals in self gravitating disks, the baroclinic instability dis- 
appears for radial wave numbers much larger than the azimuthal wave numbers {k^it) ^ ky). 
This provides a potential problem for the instability in a shear flow. The radial wave num- 
bers always tend to increase because of the winding up of spirals. In a WKB-approximation, 
one would immediately erase all baroclinic unstable modes. Thus, we can not make use of 
this simplification. 

Non-axisymmetric perturbations cannot have a simple waveform (Goldreich & Lynden- 
Bell 1965) because of the effect of the shearing background on the wave crests. Thus one 
adopts Lagrangian shearing coordinates 

r' = r, (24) 
0' = 0-fi(r)t, (25) 
t' = t, (26) 

comoving with the unperturbed shear flow. Then, it follows that, 

dr = dr'-t^d^„ (27) 

= d^,. (28) 

The new time derivative is defined as: 

dt = dt' - n{r)d^. (29) 

The azimuthal velocity is also split into two parts: 

v,p = v^,' + Q{r)r. (30) 

Linear perturbations are assumed to have the space dependence e'-iK'^' +'^<t>') Local La- 
grangian coordinates. Thus, one has now a time dependence in the radial wave number with 
respect to the non shearing one: 

Tfl 

kr ^ kr{t) ^ k'^-mt— = k'^ + q—nt, (31) 

dr r 

using D, oc r~*. One starts from the vertically integrated hydrodynamic equations in cylindri- 
cal coordinates, performs a linearisation, and solves the linear equations for the frequencies 
cu, and the growth-rate F = —icu. 
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In Lagrangian shear coordinates (dropping the primes at all coordinates), the Eqs. (3) 
to (6) transform to: 

dtT. + \[dr- t^d^^ i^^'^r) + ^d^^v^ = 0, (32) 

dtV^ + Vr {^r - ^^^'/-^ (^'Z' + + ^^'t>'"'t> + ^^"^ ^ = 

Using a local approximation simplifies the basic equations. Assuming r — R-\- (f) — y/R, 
V(j) — D,qR + Vy and x <^ R, the ideal hydrodynamical equations are in a frame corotating 
with the mean orbital frequency D,o' 

dtT. + {d^ + qn^tdy)T.v., + dy'Lvy = Q, (36) 
dtVx + Vx{dx + qVLQtdy)vx + VydyVx^- — {dx + qVtQtdy)P + 2VLQVy-2qQ.QX, (37) 

dtVy + Vx{dx + qQotdy)Vy + Vx^^O + VxdrQR + VydyVy = " ^ (9y P, (38) 

dtP + {dx + q^otdy) P + Vj;9j/P = -7P ((S^; + qilotdy) + dyVy) . (39) 

These equations are linearized, and using Eq.l3 and the definition of Oort's constant (see 
Binney & Tremaine 1987), 

P(r) = -^[a,(0r) + 0] = (|-l)0, (40) 
it follows (dropping the index at in the following) that 

dtT.' + {dx + qntdy)T.Qv', + T.odyv'y = 0, (41) 



dtv', = -^{dx + qQtdy)P' + ^dxP + 2Qv'y, (42) 
^0 ^0 



d,v'^-2Bv', = ^^dyP', (43) 
dtP' + v'^dxP = -jP{{dx + qntdy)v', + dyv'y). (44) 

The radial pressure gradient is given by d^P — —^P. This simple model helps to understand 
the basic mechanism, and its flaws will be overcome by the global model. Indeed, the global 
analysis delivers a limit on the pressure gradient /? while this local model only restricts /9 7^ 0. 
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Using the quantity B gives the possibihty to keep the analysis general for arbitrary 
rotational profiles (e.g. Keplerian rotation, rigid rotation or Rayleigh unstable fiows). 

For the perturbed state we assume: 

Q{x, y, t) = Q^e-'^^t-kyy-k.{t)x) ^ ^45^ 

even complex exponentials of the form assumed for Q{x,y,t) are not exact solutions to 
(41)-(44), and a priori one would not expect them to be good approximations unless the 
dimensionless shear rate q is small. The point is, that u generally is also a function of time 
u{t). On the other hand, following Goldreich and Lynden-Bell (1965), Balbus and Hawley 
(1992), Ryu and Goodman (1992), and many others [ultimately this technique goes back 
to Kelvin at least], one can derive from Eqs. (41)-(44) a set of four ordinary differential 
equations in time for the functions Qait) defined by: 

Q{X, y, t) = Q^{t)e'^k^^'^^^tqntky)x+ikyy _ ^4g^) 

The behavior of Qa{t) is not sinusoidal except approximately so at large \t\, where it can 
be decomposed into a sum of terms with slowly changing WKB frequencies and amplitudes. 
In this way one can calculate the correct swing amplification factor and compare with the 
approximate results in the second panel of fig. 5. This numerical check was carried out and 
is described in section 5. The results show that the approximate solution Eq. (45) leads to 
a qualitatively correct result in the analytic approach. 

Using Eq. (45) to search for eigen- values of Eqs. (41)-(44) one obtains: 

-iio-^ + i{kx + qkyVtt) v^a + ikyVya = 0, (47) 
^0 

- iujv^a - '^^Vya + i {k^ + qkyQt) ^ = 0, (48) 

Pa 

-2Bvxa - ioovya + iky-^ = 0, (49) 

^0 

/ /3\ P P Pa 

{i{k^ + qkyVtt) 7 - D + ^ky'y—Vya - ii^^ = (50) 

V -ft/ 2^0 ^0 ^0 

The two terms proportional to /3 make all the difference with respect to the stability analysis 
found in the hterature (e.g. Binney & Tremaine 1987). Without them one retains only 

the known stability criteria, because setting /? = produces a third order system out of this 
fourth order system. With (3 — surface density E and pressure P are no longer independent. 

The mean pressure can be replaced by the expression for the sound speed as: 

7^ - clix). (51) 
^0 
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+ 2r 



n {k,ky + klqnt) Cl + %ky^ (2 - I) ^ 



0. 



(52a) 



The linear modal Eqs. (47 - 50) permit solutions for F = —iu only for a vanishing determi- 
nant: 

r i {kx + qkyflt) 

M r 

— 2B J. Kvjy 

c2 (A;, + #j,l]t) - 1^) tc% r 

One obtains the dispersion relation: 



iky 
—2^2 i {kx + qkyflt) 



(53) 



If one uses the dimensionless azimuthal wave number m = kyR and the time-dependent 
dimensionless radial wave number n{t) = kxR + qmflt this relation simplifies to: 



p4 ^ ^2 



+ 2r 



n + m 



- [mn) +1771— 
2 7 



(-1) 



/3W /i/ 



7^ 



(54) 



This relation is not restricted to Keplerian disks but valid for all centrifugally supported 
rotating bodies with the general description for q. In the case of a Keplerian disk with 
q — 1.5 it is 



p4 ^ p2 



+ 2r 



- imn) -\- im , 
4 ' 74 



H 
H 

[35 



1 



/H 



R 



r 



(55) 



This 4th order equation has two high and two low frequency solutions. In Fig. 1 the 
growth-rate Re{r) and frequencies /m(r) are plotted for the case H/R — 0.1, m — 10, 
P — 1 in units of the rotational frequency Q as a function of the time-dependent radial wave 
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number n{t). The roots were found using the IDL routine FZ_ROOTS. One recognizes in 
fact two solutions with frequencies close to zero and two solutions with frequencies larger 
than the rotational frequency fl. The high frequency solutions are the ones which show 
anti-symmetric behavior in growth-rates with respect to time. Those are the ones which are 
negative for n < 0. 

The dispersion relation can be more easily interpreted by considering the high frequency 
and low frequency limits separately. 



For the high frequency part |r| ^ Q one may neglect the constant term but not the 
term linear in F. This is a depressed cubic equation and its solution is given by Cardano's 
method, which shall not be done here. For the m — case one has a 2nd order equation 



This axisymmetric solution is the usual dispersion relation of sound waves, with the epicyclic 
frequency k, and there is no growth unless < —> q > 2, which is known as Rayleigh's 
stability criterion. 



Considering the low frequency |r| <^ J7 part of the dispersion relation one neglects the 
highest order F^, thus retaining a second order system for the low frequency waves only: 



3.1. High Frequency Waves 




(56) 



3.2. Low Frequency Waves 




(57) 



which is a quadratic equation with the roots: 




WW 



- 12 - 



± — m — il 
7 R 




^2! 



(58) 



In the barotropic case (3 — there is only one solution, which is: 



10 



—m 



qii 



02 



(59) 



This is not to be confused with an instability, because it is antisymmetric with respect to 
n{t). This means there is no effective growth over a period of time: 



/to 
■to 



rio{t)dt = 0. 



(60) 



For a baroclinic (3^0 rigid rotating object {q = 0) (e.g. a star) the growth will stay 
linear for all times, because n is no longer time-dependent, i.e. there is no shear which winds 
up the spirals: 



—m- 



7 



+ rn? + 4 

/5 



± —m—il 



(61) 



This equation should be appropriate to describe baroclinic modes in centrifugally supported 
stars and other fast rotating objects. 



For a Keplerian disk {q = 3/2) one has nevertheless time-dependent radial wave numbers 



nit): 



\n{t) + i 



-m- 



5i 
4 7 



± —m—Q 
7 R 




33 ii2 



n(t)2+m2 + gg 
In the barotropic case (3 = most terms vanish and the growth-rates are 



(62) 



ro(t) 



\inn{f) 



(n(i)2 + m2) + || 



n. 



(63) 
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What is described here, is how a pattern with leading spirals gets unwound and its amplitude 
increases over this period. As soon as the spiral is open, and the winding up to trailing waves 
starts, the amplitude decreases again for the same amount it initially was increased. This 
solution has no imaginary part which means it is completely in phase with the local rotation. 
This dependence of Tq on time (respectively time dependent wave number) is shown in Fig. 
2. and Fig. 3 as the dotted lines. Fig. 2 and 3 only differ for the range of radial wave numbers 
plotted. The chosen parameters are H/R — 0.055, and m = 48 The anti-symmetry between 
negative and positive wave numbers n tells us that here effectively no amplification is taking 
place. 

In the baroclinic case /3 7^ a bifurcation of the solution branch occurs. One gets two 
different waves with growth-rates Re{T~^) and Re{T~) and with two different frequencies 
(/m(r"'"), 7m(r~), which resemble the phase velocity of the different eigenf unctions. The 
r+ solution is plotted as a sohd fine in Fig. 2 and 3 for /3 = 1. For large negative values 
of n, which corresponds to tightly wound leading waves, a behavior close to the barotropic 
solution occurs. But as the spiral wave approaches a state where it is basically open the 
amplification does not go to zero but remains small but positive for all times. During the 
up winding of trailing waves there is no decay but even a further transient amplification of 
the spiral pattern. 

The other solution branch F^, which is given by the dashed line, has zero growth as 
it starts with negative wave numbers and even decays when it approaches the unwound 
spiral n{t) = state (i.e. the respective time). As the spiral winds up to the trailing state 
this pattern has a decaying amplitude and is thus of no relevance for the global baroclinic 
instability. 

The two solutions have different eigenfunctions with different phase velocities and growth 
behavior. The important effect in the global baroclinic instability is the bifurcation of the 
solution branch into a transiently linear growing mode and a monotonically linear decaying 
mode. The idea is now that the transient linear growth can lead to the development of non- 
linear turbulence including vortices and explain the observed instability in the non-linear 
simulations of Klahr & Bodenheimer (2003). 

Hence, this shows that our vertically integrated baroclinic models have no pure linear 
unstable modes in the presence of shear. The instability appears to be transient and one has 
to investigate mode coupling or finite amplitude effects in a non linear analysis. In the case 
of rigid rotation one finds nevertheless a true linear instability. This might be of importance 
for rigid rotating objects hke stars close to their breakup velocity. 

The local analysis gives no minimum value for the entropy gradient or, correspondingly. 
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p. Any /3 7^ (positive and negative!) leads to transient growth ~ But with a global 
analysis the situation changes. 



4. Restricted Global Linear Analysis 

It is not possible to perform a general global analysis with simple k Fourier modes in 
the radial direction {k independent on r). Thus the analysis will be restricted to a special 
family of solutions which are described by 

Q(r,0,t) = i?e[Q„e-^('"*-"^'^-«"^'^('')*)]. (64) 

The assumption dtw << 1 is the same as we used in previous section. These arc spirals 
which are open at the time t — 0, leading for t < and trailing for i > 0. Sphtting the 
quantities in the global hydrodynamic equations (Eq. 32-35) in Lagrangian coordinates into 
a mean and a fluctuating part, without doing the local approximation and dropping terms 
of higher order, leads to the global hnearized equations: 

dtE' + ^(dr + ^Qtd^^{r^ov'r) + yd^v'^ = 0, (65) 

zjq V r / rzjQ 







dtP' - = -7P (^^ (dr + ^ntd^^ rvl + ^d^v'^^ , (68) 

Solution of the linearized equations can be written as a sum of terms of the form of Q' — 
Q{r')ae~^'^'^^'~"^'^'\ Without a radial variation of the initial perturbation Q{r,(j),t) at (t=0) 
(i.e. drQa — 0) one has for the radial derivative: 

Sa / 1 m \ m 

-ICO— + { - +l — qUt \ Vra + i—Vd>a = 0, (69) 

So yr r / r 

^ ;PEa-iuJVra-2nv^a+^Pa = 0, (70) 



ITTl ZTTl 
+ — qntPa-2BVra-iuJV4,a + ^Pa = 0, (71) 

P ITflP 

— {j{l + imqQt)-P)vra + 7—^v^a-i'^Pa = 0, (72) 

These four equations (Eq.(69),(70),(71), and (72)) constitute a linear system of equations 
with r = —icu and using the definition: 

n{r) = mqn{r)t (73) 
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this can be written as: 



/ r 
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/ 



(74a) 



For non trivial solutions the determinant of the coefficient matrix must vanish. Thus the 
following dispersion relation is obtained, using S = (| — l) O: 



p4 ^ p2 



+ r 



+ — in — ^ „ 

^2 



2mn^n% + 2imn ( (2-^)(^-l]% 



^4^2 



0. 



(75) 



A numerical method like FZ_ROOTS from IDL allows us to investigate the entire parameter 
space {q,n,m, P, H/R) for this general dispersion relation. The results for m = 60, g = 1.5, 
P — 1 are plotted in Fig. 4. One recognizes just hke in the local case two high and two low 
frequency solutions. The high frequency ones have an anti-symmetric behavior with respect 
to time while the low frequency ones are not at all symmetric in time. The interested reader 
is invited to explore this parameter space on his/her own. 

For an accretion disk with g = 1.5 the complete dispersion relation is: 



p4 ^ p2 



+ F 



2 I 2 • 
m + n — m 

7 

3 /5/5 
--mn + 2im — 
2 V47 



9' 



P^ 2(0^ 



0. 



(76) 



If the last term vanishes because /3 = 0, one retains a third order equation which leads to 
the Rayleigh stability criterion. This expression looks very unhandy but one can again split 
the solution into a high and a low frequency solution for F. 
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4.1. High Frequency Waves 



For the high frequency part | F | ^ f2 one neglects again the constant and the hnear term 



in F 



J- h 



-{-"n? — in — - 
+ 4 



(77) 



which is for m = basically nothing else than a Solberg-H0iland criterion for stability against 
axisymmetric perturbations: 



> gdr (InP) (V - Vad) - -^dr (r'n^) 



(78) 



For instability, j3/^ has to be steeper than (R/H)^ which is in a usual thin accretion 
disk probably never the case (e.g for a reasonable ^ ~ 0.1 one would need a P > 100 to 
make the flow unstable). Such an instability may become important if the rotational profile 
becomes steeper because then k. drops. 



4.2. Low Frequency Waves 

Considering the low frequency part of the dispersion relation one neglects Ff <S 1, thus 

(|-)(f) 

(-l)^ 



2,2- A 

m +n — m 4 



7 



+ r. 



2mn- + 2im 



T \R 

which can be solved for F; defining 



0^ = 0. 



C2 = 



2 I 2 • 

m -\-n — m 

7 



mn- + im 



^(|-)(f) 



^ _ /3W (H 
^0 ~ ~Ji 



r 



R 



(79) 



(80) 
(81) 
(82) 
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This leads to the quadratic equation 



r, = 



Ci _|_ ^/ Cf + C0C2 



n, (83) 



C2 C2 

The complete solution for Keplerian disks q — 1.5 is: 

+ — m — ^ + ^ 

Jilnr + ^|n (|f - l) - (|^ - l)V ^ (f )^ (m^ + - m - f + f ) 

± ^-^^ : ^-^2 W 

+ — — ^ + ^ 

In Fig. 5 the growth of the positive mode is plotted together with the integrated amplification 
for a disk with H/R — 0.1, m — 1, and P — 1 (solid line). The predicted amphfication of 
~ 100 fits to the numerical result (section 5) for the same parameters (see Fig. 8). This 
confirms that our simplified dispersion relation yields useful results. In addition we plotted 
the baroclinic solution for /9 = (dotted line) to show that the important effect in the 
baroclinic solution is the bifurcation of the solutions into a positive and negative solution 
branch. 

The occurrence of this bifurcation depends on (3. To study this bifurcation it is sufficient 
to focus on the behavior of the solution at n{t) = 0. 



r^C) = '" „,; (86) 



^ « ST (87) 

in ^ ^ H2 

The first part is always imaginary determining the frequencies of the waves and the second 
term leads only to a bifurcation of the modes in growth-rates, if the root has a real part. 
Thus the condition for bifurcation and by that for transient growth is 



47 J \Rj \ m 

In general the critical (5 for bifurcation is given by this third order polynomial. The depen- 
dence of r on /3 is plotted in Fig. 6 for the m = 60 modes. It appears that m ~ is 
required for maximum growth. 
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Basically there are bifurcating solutions for /3 < —0.4 which means a positive entropy 
gradient and for /3 > 0.25. The unstable regime for /3 for modes with m < R/ H is given by 

^7 < /5 < 47 (89) 
in the case of a radially constant mean density. For larger m (e.g. m > y) 

is valid. As a maximum reasonable azimuthal wave number one can assume m — 2t:R/H 
because for larger wave numbers the vertical structure can not be neglected anymore. Thus a 
critical minimum (3 for the occurrence of a global baroclinic instability should be = 0.22. 



5. Numerical Integration 

Non-axisymmetric solutions of Eqs. (41)-(44) with a time-dependence of the form exp{iujt) 
do exist only in the absence of shear {B — —Q) . To justify the usage of this simplified time 
dependence in a shearing system, as in section 3 and 4, one has to perform a numerical test 
integration of Eqs. (41)-(44) without substituting the time dependence with exp{iu!t), i.e. 

dt^a = -irb{t)v^a - imvya, (91) 
dtv,a = -pPj: + 2vya-in{t)Pa, (92) 

dtVya = --v^a- imPa, (93) 

dtPa = -{in{t)^-(3)Pv^a-im^Pvya, (94) 

with 

3 

n{t) = -mt (95) 

The integration is carried out using an implicit Runge-Kutta method (RADAU5) of order 
5 with step size control (see Hairer & Wanner 1996). For the parameters, we assume Sq = 
1.0, = 1.0, Ro = 1.0, m = 60, H/R = 0.1, and -AO/^o < t < 40/Qo to compare the 
results to that presented in Fig. 5. In Fig. 8, we plot the time evolution of the most unstable 
eigenfunction. The highest frequencies with cu » were properly resolved in the numerical 
integration, but they were removed from the all the following plots via a simple time filtering 
procedure. This does not change lower frequencies. 

The initial solution of the most unstable eigenfunction was found numerically for: 

n = -4o|m/no = -3600 (96) 
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by integrating the initial value problem Ea(0) = | — | j x 10 and -Pa(O) = ^^^(O) = 
Vya{0) = in time for the fixed n — —3600 and renormalising all amplitudes every couple of 
time steps to keep the amphtude of the density fluctuation at |Ea| = 10~^° until the solution 
converges. As one can already see from Fig. 1, only one solution shows growth for large 
negative wavenumbers. This is also the eigenfunction that shows the baroclinic unstable 
mode eventually. 



5.1. Case A /3 = 1 

In Fig. 8 we plotted the results from the numerical integration of the initial perturbation 
with (3 = 1. All values are scaled to the value of the initial surface density perturbation 
amplitude. The surface density perturbation has qualitatively the same behavior as in the 
analytic estimate (compare with fig. 5). The quantitative behavior is different in a sense that 
the swing amplification factor is about two times larger then estimated. On the other hand, 
the initial amplification of the pressure and velocity fiuctuations vanishes almost completely 
for large times. We also plot the amplitude of the potential vorticity perturbation q: 

,(0 = ^W±^-lo (97) 

Here uj{t) is the vorticity: 

Uj{t) = V X Va{t) = I %{t)Vya{t) - kyV^a{t)\ , (98) 

and must not to be confused with the frequencies of the eigenfunction. Potential vorticity 
(also called vortensity) is a conserved quantity in barotropic flows. But in baroclinic flows 
it can be generated, which can be seen in Fig. 8. The amplification in potential vorticity is 
not very dramatic compared to the initial state. Later in Case C (Fig. 10), we perform a 
calculation with vortensity generation starting from an initial state of zero local vorticity. 

Right after i = 0, high frequency waves show up in the amphtudes. These are the sound 
waves which are unfortunately also part of the solution. But their frequency diverges and 
the plotting routine filters them out successfully. 



5.2. Case B /3 = 



If one integrates the equations using the same initial values from last section but now 
with /3 = 0, the situation looks very similar for times t < (see Fig. 9). This fits perfectly 
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to the expectations from Fig. 5. In contrast to Fig. 8, now all amplitudes decay for times 
t > 0. This also agrees with the analytic estimates. 

In addition, one can now observe that vortensity is conserved, as it should be the case 
for barotropic flows. 



For the numerical integrations in Fig. 10, we did not use the numerically derived most 
unstable eigenfunction but an eigenfunction that contains initially zero vorticity (Sa(0) = 



potential vorticity is generated in the baroclinic case /? = 1 in comparison to the barotropic 
situation, where a; = at all times. This time, surface density is not amplified by as much 
as in the case of the pure eigen function (case A). 

We conclude that the usage of the simplified time dependence (~ exp{iujt)) for the 
derivations of the eigenvalues of the problem and the critical parameters for stability and 
instability is justified (see also Ryu & Goodman 1992). 



For completeness we also added a case with a negative /3, which corresponds to a positive 
pressure (entropy) gradient. As expected, the result (see Fig. 11) looks very similar to Case 
A, only the potential vorticity generation is smaller this time. Thus, we stress again: disks 
with any pressure gradient (negative and positive) can undergo this baroclinic instability. 



The linearized equations in the shearing sheet formalism allow a simple study of the 
basic instability cycle. Therefore, one neglects the radial derivative dx « dy and focuses 
on times close to t — 0, where the effect (and the bifurcation) is the clearest: 



5.3. Case C u;(i = 0) = 



X 10 and -Pa(O) = Vxa{0) = Vya{0) = 0). Thus, we can show how effectively 



5.4. Case D /3 = -0.5 



6. Explanation of the fundamental instability cycle 




(99) 



dtv'x 



(100) 



^dyP' + 2Bv', 



(101) 



'0 
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dtP' = V'^^P - ^PdyV'y. (102) 

A simple sketch will elucidate what happens in the disk (see Fig. 7). If there is an azimuthally 
confined over-density in the diskE' > 0, but no other perturbation (especially no perturbation 
of the pressure) , then equation 100 predicts a radial inward drift of matter, 

d^vJ^-^Po^ => Vx'<0, (103) 

because the disk is no longer in radial equilibrium. This is basically radial buoyancy in 
a rotating system. A region with lower density but constant pressure must have a higher 
temperature. The radial inward drift Vx < then decreases the local pressure as follows 
from Eq. 102: 

dtP' = ^Pqv' => P' < 0. (104) 
r 

But, this will result in an azimuthal velocity towards the pressure minimum (see Eq. 101): 

dtV'^ = -^dyP' => dyV'y < (105) 

which in the end will lead to a further increase of azimuthally concentrated density (see Eq. 
99): 

dt^' = -J^odxvJ => S'>0. (106) 
Now the cycle is closed, and a further amplification can follow. 



7. Discussion and Conclusions 

7.1. Do Accretion disks have a negative Entropy gradient P > 0.22 ? 

Isothermal disks have no entropy gradient and will be stable. But, not every disk can be 
assumed to be isothermal. Three scenarios in which the instability criterion can be fulfilled 
are given here: 

• A Disk close to the central object or to the boundary layer: 

In such a situation, the temperature gradient is determined by the radial radiation 
transport in the disk. The temperature profile can be very steep in an optical thick 
disk, irradiated only at the inner edge. 

• A protostellar disk: 

A disk that has just formed after the collapse of a molecular cloud core is also likely 
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to have a temperature gradient slightly steeper than 0.5 (see Fig. 4. in Yorke & 
Bodenheimer 1999). 

• A disk with self sustained accretion: 

Under the assumption of a constant accretion rate, the locally produced energy gives 
a temperature profile of T ~ which will be unstable to the Global Barochnic 

Instability. Two- and three-dimensional radiation hydro simulations of accretion disks 
showed a temperature slope of T ~ while the surface density was constant over 
wide ranges of radii at 3.5 — 6.5 AU (Klahr & Bodenheimer 2003). 

Detailed studies on the actual entropy gradient in various disks will have to follow this study. 

So far only the simple case of constant surface density was studied in the linear analysis, 
which is indeed the proper solution for an a-model with dust opacities for a constant accretion 
rate at about 5 ± 3AU (see Bell, Cassen, Klahr, & Henning 1997). 

Nevertheless, an analysis for general f]^, will have to follow. What one can expect 
from this more general analysis is a stability criterion that not only depends on the entropy 
gradient but the entropy gradient in relation to the local pressure gradient in the fashion of 
the Schwarzschild-Ledoux criterion: 

dP dS ^ (107) 
dr dr 

Only when entropy gradient and pressure gradient point in the same direction a convective 
instability occurs. As long as the surface density is radially constant, this criterion is always 
fuimUed. 



7.2. Closing remcirks 

Rotating centrifugally supported systems (e.g., accretion disks, galaxies, fast rotating 
stars) can be hydrodynamically unstable to non-axisymmetric perturbations under the in- 
fluence of a radial entropy gradient. Even if they are stable with respect to the Rayleigh and 
S0lberg-Hoiland criterion, the radial entropy gradient can produce a Baroclinic instability 
forming slowly oscillating waves. A similar effect was known in the case of planetary surfaces 
(Classical Baroclinic Instability) and has now been demonstrated for accretion disks. 

While the instability produces linear growth for rigidly rotating bodies, which is known 
since the 1940s (e.g. Cowhng 1951), it is weakened by the radial shear inherent in accretion 
disks. Thus, in sheared geometries, the instability is transient, leading only to a finite 
amplification in the linear regime. Non-linear effects can then possibly pick up and lead to 
further amplification. 
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It could be possible that the nonlinear vortices arise through the following three steps. 
First, swing amplification of an initially nonvortical leading perturbation produces a vortical 
trailing one. We know that in a barotropic disk, a perturbation in vorticity per unit surface 
density is permanent (absent dissipation) once it is created. In these linearized calculations, 
one expects that the trailing vortical perturbations will be approximately constant in ampli- 
tude at late times because the baroclinic source terms for vortensity become small when the 
wave is tightly wrapped (large kx{t)/ky). There will then be an almost axisymmetric pattern 
of radial maxima and minima in vortensity, which gradually become more pronounced as 
kx{t) increases. Thus, the second step may be a Kelvin- Helmholtz instability of this pattern, 
resulting in small-scale vortices ("Kelvin's cat's-eyes"). The third and final step would be 
the merging of small anticyclonic vortices into larger ones. 

This would explain the results from non-linear simulations (Klahr & Bodenheimer 2003) 
and also give a physical model for the vortex formation in protoplanetary disks, which is 
supposed to be the starting point of planet formation (Klahr 2003). It would be useful to 
make predictions for the shape, size, and frequency of these large-scale anti-cyclonic vortices 
before they become observable with interferometric instruments like ALMA (Wolf & Klahr 
2002). 
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Fig. 1. — Solutions of the local dispersion relation: Growth rates and frequencies of the 
m — 10 modes as a function of the time-dependent radial wave number n{t) = no + 3/2ftm.t 
in a narrow range from —20 < n < 20. For explanation see text. 
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Fig. 2. — Solutions of the local dispersion relation for /3 = 1: Growth rates and frequencies of 
the m = 48 modes as a function of the timedcpcndent radial wave number n{t) = no+3/2Qmt 
in a narrow range from —10 < n < 10. The dotted line is the barotropic case Tq, while the 
baroclinic case is given by the solid r+ and dashed T~ curve. For explanation see text. 



-27- 



[rowth 



1.0 



0.5 



0.0 



-0.5 



-1.0 

-1000 -500 500 1000 
n(t) 



frequency 



111,111 


, 1 1 1 1 , , 1 1 




[ 




/' 




' / 




! / 




\j 







0.000 
-0.005 

-0.010 

-0.015 

-0.020 

-0.025 



0.030 

-1000 -500 









1 / 




/ 

/ : 

1 : 

1 

1 




1 : 
1 

1 - 

1 

1 

1 ■ 



500 1000 



n(t) 



Fig. 3. — Solutions of the local dispersion relation for (3 = 1: Growth rates and frequencies of 
the m = 48 modes as a function of the timedependent radial wave number n{t) = nQ+3/2flmt 
in a wide range from —1000 < n < 1000. The dotted hne is the barotropic case Fq, while 
the baroclinic case is given by the solid F+ and dashed F~ curve. Same data as Fig. 2. For 
explanation see text. 
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Fig. 4. — Solutions of the global dispersion relation for (3=1: Growth rates and frequencies 
of the m = 60 modes as a function of the time-dependent radial wave number n{t) = 3/2fl'mt 
in the range from —40 < n < 40. 
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Fig. 5. — The transient unstable mode of the global dispersion relation: Growth rates F"*" 
and integrated amplification Q{t) of the m = 60 mode as a function of time in units of orbital 
periods. Compare the baroclinic case (/3 = 1: solid line) with the barotropic case Fq {(3 ~ 0: 
dotted line). (Unlike in the previous figures time t is used instead of n{t) — 3/20,mt for this 
plot.) 
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Fig. 6. — Solutions of the global dispersion relation: Growth rates and frequencies of the 
m = 60 modes at time t = as a function of the radial pressure gradient f3. Bifurcation in 
growth-rates only occurs for certain values of (3. The dotted Une is the stable m = case, 
while the barochnic case is given by the sohd r+ and dashed T" curve. For explanation see 
text. 




Fig. 7. — Schematic process of the Instabihty cycle. For explanation see text. 




Fig. 8. — Evolution of the perturbations in a Keplerian disk for Case A (f3 — 1): The 
amplitudes of surface density S^, pressure Pa. root mean square velocity Vrms and potential 
vorticity q are plotted. For explanation see text. 



-33- 




Fig. 9. — Evolution of the perturbations in a Keplerian disk for Case B (/3 = 0): The 
amplitudes of surface density S^, pressure Pa. root mean square velocity Vrms and potential 
vorticity q are plotted. For explanation see text. 




Fig. 10. — Evolution of the perturbations in a Keplerian disk for Case C {P — 1 and zero 
initial vorticity): The amplitudes of surface density Ea and potential vorticity q are plotted. 
For explanation see text. 
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Fig. 11. — Evolution of the perturbations in a Keplerian disk for Case D {P — —0.5): The 
amplitudes of surface density Sa, pressure Pa. root mean square velocity Vrms and potetntial 
vorticity q are plotted. For explanation see text. 



